Dynamic patterns of electroosmosis peristaltic flow of a Bingham fluid model in a complex wavy microchannel

The purpose of this paper is to present a rigorous analysis of streamline patterns and their bifurcation to a viscoplastic Bingham fluid model that involves heat and mass transfer in an electroosmotic flow through a complex wavy microchannel. The Bingham fluid act as a solid medium in the core layer, which divides the channel into three distinct sections utilized to model the problem as a switched dynamical system between these zones. To track multiple steady states (stagnation points) and related trapping phenomena, we perform both analytical and numerical bifurcation analysis of each subsystem with respect to different physical effects such as electrical double layer thickness and Helmholtz-Smoluchowski velocity. The key feature of the technique presented here is its ability to reveal the peristaltic transport characteristics of the Bingham fluid model in the presence or absence of symmetric flow properties. The primary novelty here is the ability to regulate the location and stability of the equilibrium points in the domain of interest. This leads to the detection of global bifurcations that reflect important dynamic elements of the model. Our results highlighted a new category of complex behavior that controls transitions between qualitatively different transport mechanisms, as well as a class of non-classical trapping phenomena.

www.nature.com/scientificreports/ Fig. 1. Viscous dissipation and joule heating are also considered. A mathematical model of wall deformation 17,18 is shown below.
where ε i needs to fulfil the condition m i=1 ε i ≤á. Bingham fluid is a viscoplastic fluid, a type of non-Newtonian fluid in which the flow field is divided into two regions: an un-yielded zone in which the fluid is at rest or undergoes stiff motion, and a yielded zone in which the fluid flows like a viscous liquid. In the un-yielded zone, the second invariant of the extra stress tensor is less than or equal to the yield stress and a constitutive relation is undefined. In the yielded region, this invariant exceeds the yield stress and a constitutive relation exists for the extra stress tensor. Thus, the location and shape of the yield surface(s), i.e. the interface between these two sets, is also a part of the solution of flow problems of such fluids. Viscoplastic fluids occur in various chemical, metal, and food industries, e.g., margarine, mayonnaise and ketchup. The constitutive equation of an incompressible Bingham fluid is based on the assumption that the fluid remains at rest or moves as a rigid body if the second invariant of the extra stress tensor S xy is less than or equal to the yield stress S 0 .
The constitutive equation 18,25 that is required to describe a Bingham model is: where S 0 represents yield stress and γ represents the rate of strain tensor. As a direct consequence, the Bingham fluids act as a solid medium in the core layer.
Bingham fluid applications. The hypothetical viscous fluid has a yield strength that needs to be exceeded before it may flow. A lava channel, for example, is defined as a stream of flowing lava contained within zones of static (i.e., solid and motionless) lava or lava levees. Levees may not exist in the original channel until the parental flow settles over what forms the channel and produces natural levees. Therefore, lava behaves as a multiphase non-Newtonian fluid, for example see 40 . In this context, most lava flows are obvious applications of Bingham fluids. additionally, Bingham Fluid can be designed as the interaction of blood's non-Newtonian nature and its flow through arteries, such as microvascular blood flow through a complex wavy microchannel, see 10,18 .
Governing equations. The governing equations that serve as guiding principles for electroosmotic fluid flow 18 are given by: (2) is the differential operator. The variables and parameters used in the preceding system are defined in List of symbols section. Due to the presence of an electric double layer (EDL) in the microchannel, the electric potential is calculated using the Poisson equation 12 , which can be expressed as: such that ρ e = ze n + − n − .
In this context, the Nernst-Planck equation for ionic number distribution is utilized to assess potential distribution as follows: We introduce the following set of non-dimensional variables and parameters to help get explicit analytical solutions to the governing equations.

Explicit analytical solutions
Debye-Huckel modifies the Poisson equation as follows: where, D l = aze 2n 0 εT a k B . Additionally, the dimensionless quantities (10) allow for a reduction of the Nernst-Planck Eq. (9) to establish the ionic distribution as follows: Solving the above equation with the associated conditions ( ∂n ± ∂y = 0 at ∂� ∂y = 0 ), and ( n ± = 1 at = 0 ) we obtain: Thus, the Eq. (11) becomes (10) www.nature.com/scientificreports/ The above equation can be linearized using low-zeta potential approximation (i.e, sinh (�) ≃ � ), as: The potential function is obtained as a explicit solution of Eq. (15) subject to the boundary conditions ( 0 = ∂� ∂y | y=0 and 1 = | y=h ) as: The governing Eqs. (1)-(7) for electroosmotic fluid flow are reduced to non-dimensional forms using the dimensionless quantities (10), lubrication theory, low Reynolds, and large wavelength approximations, as follows: Physical boundary conditions for temperature, concentration, and velocity are imposed as 18 : The axial velocity is determined by solving the Eq. (20) with boundary conditions (24) as follows: The solution for the normalized temperature distribution is found by solving Eq. (22) and using the associated boundary conditions (24) as follows: where (to simplify, use D p = ∂p ∂x ), (14) ∂ 2 � ∂y 2 = D l 2 sinh (�). where α 1 = k r S c and α 2 = k r S c N d .

Bifurcations of equilibria and nonlinear behavior
The main goal of this section is to describe and control the nonlinear behavior of all flow modes using dynamical system theory and state space simulations. The equations of motion of individual fluid particles in a two-dimensional flow can be written in classical form, with the stream function acting as the Hamiltonian switched dynamical system. The following regions are defined based on the problem formulation and stage configuration flow.
(a) Upper region of non-plug flow 1 : Thus, for the vector ξ ∈ R 2 , the switching lines are defined by y = h pl and y = 0 and the vector fields f i (ξ , ϑ) are smooth functions on the corresponding region i , i=1,2,3. The switching system is characterized as: It should be emphasized that no dynamics occur for any ξ ∈ � 2 since the fluid velocity is considered to be constant (i.e., f 2 (ξ , ϑ) = c ) throughout any cross-section of the channel perpendicular to the channel axis, which means that all particles in a given cross-section in 2 have same velocity and direction of motion. Further, for all ξ ∈ � 3 the vector filed The equilibrium points (or stagnation points) are the solution points ξ in a flow field where the local velocity ξ of the fluid found to be zero. The classification of equilibria into admissible, virtual, and boundary points is critical for capturing the system's dynamic.
Definition 1 Assume that P ∈ R 2 is a equilibrium point of the system (27), then (a) If P ∈ i and F j |P= 0 for any i = j , i, j = 1, 2, 3 , then P is referred as an admissible(valid) point. (b) If P ∈ i and F j |P= 0 for any i = j , i, j = 1, 2, 3 , then P is referred to as a virtual point (because it is not located in its associated region). (c) If P := {x ∈ R 2 | y = h pl and f 1 (P) = 0} or P := {x ∈ R 2 | y = 0 and f 3 (P) = 0} , then P is referred to a boundary point.
BrD p , www.nature.com/scientificreports/ Switches between distinct types of dynamical behavior may be performed by describing the transition of virtual to admissible equilibrium points, as well as the possibility of bifurcation in each domain independently.
It is beneficial to consider a specific situation in order to fully comprehend the basic flow motion through the wavy channel. Hence, we consider the flow mode when ú e = 0 and therefore the two-dimensional flow is defined by stream function as: Note that if m ≥ 3 , it is difficult to find the zeros of the function (31) explicitly, we can compute the values of x using various numerical methods, such as the Newton method. Figure 2 shows that the nonlinear behavior of the function (31), whose roots are observed as x-intercepts, occurs at x-values where the function value is 0.
The configuration of the possible bifurcation of the above explicit formulas of equilibrium point is sufficient to report the most of the significant different flow behavior that may occur in this situation. It is clear in Fig. 3b that the center points generate fluid boluses around them, and the saddle points connection in the context of the existing heteroclinic curve forms a trapping zone.
In general situation, ú e � = 0,  We are unable to proceed with the explicit bifurcation analysis of the above vector fields due to their high level of nonlinearity. Therefore, in the following section, we will extend our analysis of bifurcation using an efficient numerical technique.

Numerical bifurcation analysis. Numerical analysis is an important tool in dealing with nonlinear bifur-
cation problems in various biological and physical systems. One of the main advantage of this approach is that it is used to measure the steady-state curves for an undetermined system of equations.
The solutions to the following piecewise nonlinear system are helpful in determining all potentially admissible, virtual, boundary, and degenerate equilibrium points, as well as their bifurcation.
A nonlinear system F x, y, ϑ = 0 may have an infinite or finite number of roots, and these roots can be extremely sensitive to small changes in one or more parameters of ϑ . Therefore, the main idea is to present an optimal technique for computing all roots of F x, y, ϑ and detecting various types of behavior around these roots, resulting in a complete picture of various multiphase (i.e. in 1 and 3 ) flow behaviors.
It should be noted that there is no equilibrium point for the system ξ = f 2 x, y, ϑ , y ∈ � 2 . However, some virtual points for the system (35) can be located in the domain of 2 (see for example Fig. 3b). As a first stage in computing the equilibrium points, a fixed is chosen to be a sector-like domain, which is defined as: This implies that the sector shears the same domain along the x-axis, weher a, b ∈ Z . According to numerical domain decomposition, limiting the size of the solution domains improves computational efficiency and reduces the amount of computing effort necessary to solve the system (35). Then these real intervals of interest 1 and 2 will further divide into finite sub-intervals based on nonlinear system behaves at different regions (i. e., the sub-intervals are not necessarily equidistant) as: 1 = [ m−1 1 , m 1 ] and 2 = n−1 2 , n 2 , m, n ∈ Z + . Then, in each subinterval 1 and 2 , we use Newton's iterative technique to solve the two system (35) independently, and the explicit equilibrium points (30) are considered an initial guess. Finally, we use the Matlab package (Mat-Cont)for numerical bifurcation analysis to classify the output equilibrium points and calculate their linear stability and related bifurcation 39 ).

Results and discussion
In this section, we discuss the impact of different physical parameters on fluid flow behaviors using numerical bifurcation analysis.
Trapping phenomenon: Based on a dynamical examination of flow behavior around equilibrium points, this phenomenon can be discovered in two scenarios, namely local and global modes of system behavior. The flow is trapped in terms of local mode when the system (35) has center points (closed orbits) that cause the streamline to split to enclose a bolus of fluid particles. In general, the global mode is defined as structural changes in the characteristic frame that cannot be identified when examining the stability of equirbira by computing the eigenvalues of the associated Jacobian. A heteroclinic connection is a type of global mode formed when a streamline connects two or more saddle-points. One advantage of using heteroclinic connections in such a system is that the flow is completely trapped between the connection orbits. Furthermore, the distance between two saddle points for a heteroclinic connection is used to determine the minimum and maximum limits of the trapping zone. In all subsequent computations, we consider the amplitude of the various waves to be ε = � 5 j=1 3.4j 100 . In this context, Fig. 4 depicts the largest invariant curve connecting two distances of saddle points, which establishes the maximum trapping zone. www.nature.com/scientificreports/ The plug flow region 2 is identified and controlled by the parameter h pl . Hence, we discuss how changes in the parameter h pl affect the flow behavior in the upper region of the channel due to changes in the location and stability of the equilibrium points. At h pl = 0.1 and keeping all the parameters fixed, as shown in Fig. 4. In Fig. 5a, we notice that the flow behavior in region 1 changes due to the presence of a plug flow region and the heterocilinic connection between the upper and lower regions of the channel is destroyed, resulting in two distinct trapping regions, although the flow behavior in region 3 remains unchanged. When the parameter h pl is increased to h pl = 0.3 , two distinct heterocilinic connections emerge in the upper region, and the size of the trapping bolus is reduced, see Fig. 5b.
By taking the Debye length D l (i.e., the electrical double layer thickness) as a bifurcation parameter (all other parameters fixed as ú e = 1, h pl = 0.2, q = −0.18366 ), It is seen in Fig. 6 increasing the value of D l causes some equilibrium points to vanish (disappear), while others points move to become boundary points. The fixed parameters of the system are given as: ú e = 1, h pl = 0.2, q = −0.18366 , where Fig. 6a is created with D l = 3.0 while Fig. 6b is created with D l = 5.0 . This explains why the boluses are initially reduced in size and then vanish as the value of D l increases, i.e., the trapping zone does not occur with a large value of D l .   Fig. 7a and b, but the configuration of bolus dynamics inside this zone differs depending on whether ú e is positive or negative. When ú e = −5 , the size of boluses increases noticeably, whereas when ú e = 5 , the size of boluses decreases. When h pl = 0.2 , the flow in channel becomes asymmetric on both sides of the x-axis, see Fig. 7c and d. Further, the upper row of boluses converts to a family of heteroclinic connections that form a family of trapping zones that vary in size depending on whether ú e is negative or positive, whereas the lower row of boluses does not change.
Temperature fields in 1 and 3 are determined of an interesting situation when the flow is trapped in both regions using the same parameters as shown in Fig.5 (i.e., ú e = D l = 1, h pl = 0.2, q = −0.18366 ). The contour of the temperature field across the microchannel is shown in Fig. 8a-c for negative, zero, and positive values of Joule heating term Ǵ and the Brinkman number is fixed Br = 1.0 ( whereas Fig. 8d holds for Ǵ = 1.0 and Br = 0.5 ) as the clarification of the changing and periodic temperature fields evolves to a nation of periodic dynamics. Fig. 9 depicts the fluid concentration behavior when the Schmidt number, chemical reaction, and concentration difference parameters are varied.

Conclusions
The qualitative aspects of an electroosmotic peristaltic Bingham fluid model, specifically geometrical properties of flow fields such as bifurcation and the stability of its streamline patterns, are investigated. This model simulates the effect of heat and mass transfer on Bingham fluid flow through a complex wavy microchannel influenced by electroosmosis. The following are the key results of the current study: • Analytical and numerical bifurcation analysis is used to systematically identify dynamic behavior and characterize fluid flow to reveal associated physical phenomena. • Our results indicate that heteroclinic connections to saddle points are the primary cause of the trapping phenomenon in two scenarios, namely the presence or absence of symmetric flow. • The non-uniform geometry caused by the plug region and varying amplitude ratio parameters has a significant impact on the trapping phenomenon. • We assert that the trapping zone does not exist with a large Debye length D l because the model's equilibrium points vanish (disappear) or become boundary points. • The impact of Helmholtz-Smoluchowski velocity ú e has been demonstrated in a case with interesting and complex behavior. The flow, for example, generated a maximum trapping zone; our bifurcation findings explain why the configuration of the bolus dynamics within this zone varies depending on whetherú e is positive or negative.